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Feature Lines for Illustrating Medical Surface 
Models: Mathematical Background and Survey 


Kai Lawonn and Bernhard Preim 


Abstract This paper provides a tutorial and survey for a specific kind of illustrative 
visualization technique: feature lines. We examine different feature line methods. 
For this, we provide the differential geometry behind these concepts and adapt this 
mathematical field to the discrete differential geometry. All discrete differential ge¬ 
ometry terms are explained for triangulated surface meshes. These utilities serve 
as basis for the feature line methods. We provide the reader with all knowledge to 
re-implement every feature line method. Furthermore, we summarize the methods 
and suggest a guideline for which kind of surface which feature line algorithm is 
best suited. Our work is motivated by, but not restricted to, medical and biological 
surface models. 


1 Introduction 


The application of illustrative visualization has increased in recent years.The princi¬ 
ple goal behind the concept of illustrative visualization is a meaningful, expressive, 
and simplified depiction of a problem, a scene or a situation. As an example, running 
people are represented running stickmans, which can be seen in the Olympic games, 
and other objects become simplified line drawings, see Figure[T] More complex ex¬ 
amples can be found in medical atlases. Most anatomical structures are painted and 
illustrated with pencils and pens. Gray’s anatomy is one of the famous textbooks for 
medical teaching. Most other textbooks in this area orient to depict anatomy with 
art drawing, too. 
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Other than simplified representa¬ 
tion, illustrative visualization is not 
restricted to these fields. Illustrative 
techniques are essential for focus- 
and-context visualizations. Consider a 
scene with anatomical structures and 
one specific (important) structure. The 
specific structure may be strongly re¬ 
lated to the surrounding objects. There¬ 
fore, hiding the other objects is not a 
viable option. In contrast, depicting all 
structures leads to visual clutter and optical distraction of the most important struc¬ 
tures. Focus-and-context visualization is characterized by a few local regions that 
are displayed in detail and with emphasis techniques, such as a saturated color. Sur¬ 
rounding contextual objects are displayed in a less prominent manner to avoid dis¬ 
traction from focal regions. Medical examples are vessels with interior blood flow, 
livers with inner structures including vascular trees and possible tumors, proteins 
with surface representation and interior ribbon visualization. Focus-and-context vi¬ 
sualization is not restricted to medical data. An example is the vehicle body and the 
interior devices. The user or engineer needs the opportunity to illustrate all devices 
in the same context. 

There are numerous methods for different illustration techniques. This survey is 
focused on a specific illustrative visualization category: feature lines. Feature lines 
are a special group of line drawing techniques. Another class of line drawing meth¬ 
ods is hatching. Hatching tries to convey the shape by drawing a bunch of lines. 
Here, the spatial impression of the surface is even more improved. Several methods 
exist to hatch the surface mesh, see lfT5ll20ll28llJ7ll5ni5^ . In contrast, feature lines 
try to generate lines at salient regions only. Not only for illustrative visualization, 
feature lines can also be used for rigid registrations of anatomical surfaces l44l or 
for image and data analysis in medical applications IfT^ . The goal of this survey is 
to convey the reader to the different feature line methods and offer a tutorial with all 
the knowledge to be able to implement each of the methods. 








(a) 


(b) 


(c) 


(d) 


Fig. 1 Visual abstraction of the four Olympic 
disciplines: archery, basketball, football and 
handball in the style of the pictograms of the 
Olympic Games 2012 in London. 


Organization. 

We first give an overview of the mathematical background. In Section we in¬ 
troduce the necessary fundamentals of differential geometry. Afterwards, we adapt 
these fundamentals to triangulated surface meshes in Section]^ Section |^discusses 
general aspects and requirements for feature lines. Next, we present different feature 
line methods in Section]^ and compare them in Section]^ Finally, Section]^ holds 
the conclusion of this survey. 
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Fig. 2 The basic elements for differential geometry. A parametric surface is given and the partial 
derivatives create the tangent space. 


2 Differential Geometry Background 

This section presents the fundamentals of differential geometry for feature line gen¬ 
eration, which will be crucial for the further sections. We present the basic terms 
and properties. This section is inspired by differential geometry books Il5ll6l l2^ . 


2.1 Basic Prerequisites 


A surface /; / c ^ is called a parametric surface if / is an immersion. An im¬ 
mersion means that all partial derivatives are injective at each point. The further 
calculations are mostly based on the tangent space of a surface. The tangent space 
Tpf of / is defined as the linear combination of the partial derivatives of /; 

r,/:=sp,„(f| ,f| 1. 

'OXilx=u 0X2'"^=^’ 


Here, span is the space of all linear combinations. Formally: span{vi,V 2 ) {crvi + 
j3\2 1 a.P e R). With the tangent space, we can define a normalized normal vector n. 
The (normalized) normal vector n(u) at p = /(u) is defined such that for all elements 
V € Tpf the equation (v,n(u)) = 0 holds, where (.,.) denotes the canonical Euclidean 
dot product. Therefore, n(u) is defined as: 


df\ 

dx\ \ 

x^l 

lx=u dX 2 \ 

x=u 

II 


II 

II 1 

lx=u ^->^2 1 

x=u’* 
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This map is also called the Gauss map. Figure [^depicts the domain of a parametric 
surface as well as the tangent space Tpf and the normal n. 


2.2 Curvature 

The curvature is a fundamental property to identify salient regions of a surface that 
should be conveyed by feature lines. Colloquially spoken, it is a measure of how far 
the surface bends at a certain point. If we consider ourselves to stand on a sphere 
at a specific point, it does not matter in which direction we go, the bending will be 
the same. If we imagine we stand on a plane at a specific point, we can go in any 
direction, there will be no bending. Without knowing any measure of the curvature, 
we can state that a plane has zero curvature and that a sphere with a small radius 
has a higher curvature than a sphere with a higher radius. This is due to the fact that 
a sphere with an increasing radius becomes locally more a plane. Intuitively, the 
curvature depends also on the direction in which we decide to go. On a cylinder, we 
have a bending in one direction but not in the other. Painting the trace of a walk on 
the surface and view it in 3D space, we could treat this as a 3D curve. The definition 
of the curvature of a curve may be adapted to the curvature of a surface. The adaption 
of this concepts is explained in the following. Let c: / c R ^ be a 3D parametric 
curve with ||^|| = 1 . The property of constant length of the derivative is called arc 
length or natural parametrization. One can show that such a parametrization exists 
for each continuous, differentiable curve that is an immersion. So, if we want to 
measure the size of bending, we can use the norm of the second derivative of the 
curve. Therefore, the (absolute) curvature K(t) at a time point t is defined as: 

K{t) = 

To determine the curvature on a cer¬ 
tain point of the surface in a specific 
direction, we can employ a curve and 
calculate its curvature. This approach 
is imperfect because curves that lie in 
a plane can have non-vanishing curva¬ 
ture, e.g., a circle, whereas we claimed 
to have zero curvature on a planar sur¬ 
face. Therefore, we have to distinguish 
which part of the second derivative of 
the curve contributes to the tangent 
space and which contributes to the nor¬ 
mal part of the surface. Decomposing 
the second derivative of the curve into 
tangential and normal part of the sur¬ 
face yields: 


|c"(f)||. 



Fig. 3 The curve’s second derivative is decom¬ 
posed into the tangential and normal part. 
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c"(t) = projy-^y c"(f) + <c"(f),n)n , 

tangential part normal part 

where c(f) = p and prqj^x means the projection of the point x onto the space E, see 
Figure The curvature Xc(p) of the surface at p along the curve c is defined as the 
coefficient of the normal part; 


Kcip) = <c"(f),n>. 


( 1 ) 


Hence, we know that c'{t) 6 Tpf and (c'(f),n) = 0. Deriving the last equation 
yields; 


We obtain 


-<c'(0,ii) = 0 
at 

^<c'(f),n) = <c"(f),n> + <c'(f),^>. 
at ot 




Combining this equation with Equation yields 


( 2 ) 


Thus, the curvature of a surface at a specific point in a certain direction can be 
calculated by a theorem by Meusnier. We call the vectors v,w at p the maximal/ 
minimal principle curvature directions of the maximal and minimal curvature, if 
Xy(p) > Xyz/p), Xw(p) < Xv/(p) for all directions v' e Tpf. If such a minimum and 


maximum exists, then v and w are perpendicular, see Section 2.5 for a proof. If we 
want to determine the curvature in direction u, we first need to normalize u, v, w and 
can then determine a:u(p) by; 


/fu/p) = (u, v/^Xyjp) + <u,w)2x„(p). 


( 3 ) 


The coefficients of the curvature are the decomposition of the principle curvature 
directions with the vector u. 


2.3 Covariant Derivative 

The essence of the feature line generation is the analysis of local variations in a 
specific direction, i.e., the covariant derivative. Therefore, the covariant derivative is 
a crucial concept for feature line methods. We consider a scalar field on a parametric 
surface ip ; /(/) —> R. One can imagine this scalar field as a heat or pressure (as well 
as a curvature) distribution. The directional derivative of ip in direction v can be 
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Fig. 4 Given: a scalar field in the domain. Determining the gradient and using it as coefficient for 
the basis tangent vectors leads to the wrong result (grey). Balancing the distortion with the inverse 
of the metric tensor yields the correct gradient on the surface (black). 


written as Dy(p and is defined by: 

Dy(p{x) - lim 
/«—>0 


(p{x + h\) - ifi(x) 
h 


If (f is differentiable at x, the directional derivative can be simplified: 

Dyifix) = <V(^(x),v), 


where V denotes the gradient. The gradient is an operator applied to a scalar held 
and results in a vector held. When we want to extend the dehnition of the deriva¬ 
tive to an arbitrary surface, we hrst need to dehne the gradient of surfaces. In the 
following, we make use of the covariant derivative. The standard directional deriva¬ 
tive results in a vector which lies somewhere in the 3D space, whereas the covariant 
derivative is restricted to stay in the tangent space of the surface. The gradient is a 
two-dimensional vector. Actually, we need a three-dimensional vector in the tangent 
space of the surface. Here, we employ the gradient and use it as coefficients of the 
tangential basis. Unfortunately, this leads to wrong results because of the distortions 
of the basis of the tangent space, see Figure]^ The basis is not necessarily an orthog¬ 
onal normalized basis as in the domain space and can therefore lead to distortions 
of the gradient on the surface. 

One way to calculate this vector is to use the plain scalar field ip: —> R. Af¬ 
terwards, we are able to attain the gradient in three-dimensional space and project it 
on the tangent space. However, we want to use the gradient oip: R^ ^ R in the do¬ 
main of a parametric surface and compensate the length distortion such that we can 
use it as coordinates with the basis in the tangent space. One important fact is when 
multiplying the gradient with the /-th basis vector, one obtains the partial derivative 
of (p with Xi. Hence, we know that the three-dimensional gradient V(p lies in the tan¬ 
gent space. Therefore, it can be represented as a linear combination of ^ with 
coefficients a,/?: 
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V^£) = a- —+^-—. 

OXI 0X2 

Multiplying both sides with the basis vectors and using the relation ^ = (V^, 


we obtain an equation system: 


’ dxi > 


' dp ' 

/ 

dx\ 


dp 


^8x2' 

V 


“ ^dxi’ dx, ' ^dx, ’ dx2' 

ry.(K K\ + R.(K £L\ 
“ ^dxi’ dx2' P ^dx2’dx2 0 


,oj_ 0]_. ,0]_ aj ' 

' dx] ’ dx\ ' ^ dx\ ’ dx2 ' 

, df df . ,df df . 

y^dxi ’ dx2^ ^dx2' dx2^) 


g-= 


The matrix g is called the metric tensor. This tensor describes the length and area 
distortion from to the surface. The last equation yields the coefficients a,p when 
multiplied with the inverse of g: 



-1 


dxi 

dip 

V 3x2 ^ 


This leads to a general expression of the gradient for a scalar field ip: R” —> R: 



dxi ’ 


( 4 ) 


where g'^ is the /, 7 -th matrix entry from the inverse of g and ^ means the basis. 
Now, we are able to determine the covariant derivative of a scalar field by first 
determining its gradient and afterwards using the dot product: 


D„(p - (V^,w). 


2.4 Laplace-Beltrami Operator 


The Laplace-Beltrami operator is needed for a specific feature line method and will 
therefore be introduced. The Laplace operator is defined as a composition of the 
gradient and the divergence. When interpreting the vector field as a flow field, the 
divergence is a measure of how much more flow leaves a specific region than flow 
enters. In the Euclidean space, the divergence div of a vector field 0: R” —> R" is 
the sum of the partial derivatives of the components 0 ,: 


dlV0= > 


The computation of the divergence for a vector field 0: R" ^ R” in Euclidean space 
is straightforward. However, for computing the divergence to an arbitrary surface 
we have to be aware of the length and area distortions. Without giving a derivation 
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of the divergence, the components cf>, of the vector field have to be weighted by 
the square root of the determinant -y/j^ of the metric tensor g before taking the 
derivative. The square root of the determinant of g describes the distortion change 
from the Euclidean space to the surface. Formally, the divergence of a vector field 
0: R” ^ R" with a given metric tensor g is given by: 


div0 = 




( 5 ) 


Given the definition of the gradient and the divergence, we can compose both op¬ 
erators to obtain the Laplace-Beltrami operator Jgp of a scalar field (p: R" ^ R on 
surfaces: 


J<p - divV<p - 



( 6 ) 


2.5 Shape Operator 


In Section 2.2 we noticed that the curvature of a parametric surface at a specific 


point p in a certain direction can be determined by Equation]^ 


Actually, this means that the curvature in the direction c'{t) is a measure of how 
much the normal changes in this direction, too. Given is v 6 Tpf with p = /(u) and 
|v| = 1. Then, we determine the coefficients a,p of v with the basis 



-1 


(/y ^)'\ 

(v 

6x2^ 


We use (a,jS) to determine the derivative of n along v by using the two-dimensional 
curve c(t) = u -i- and calculate: 


d 

D^n := —n(c(0). 
dt 

We define S (v) := -Dyii. This linear operator is called Shape Operator (also Wein- 
garten Map or Second Fundamental Tensor). One can see that S(^) — holds. 
Note that this operator can directly operate on the 3D space with a three-dimensional 
vector in the tangent space, as well as the 2D space with the coefficients of the basis. 
Therefore, it can be represented by a matrix S. Recall Equation]^ we substitute c' 
with v and ^ by 5 v: 
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%(p) = <v,5v). 

We want to show that the principle curvature directions are the eigenvectors of S. 
Assuming vi, V 2 e R.^ are the normalized eigenvectors with the eigenvalues Ai > /I 2 . 
Every normalized vector w can be written as a linear combination of vi,V 2 : w = 
avi +/3 v2 with ||w|| = Havi +;0V2|| = +0^ + 2af5{\i,\2) - 1- Therefore, we obtain; 

'^w(p) = <w,5w> = ^[(a^-yS^)(Ti -/t 2 ) + Ti +T 2 ]- (7) 

One can see from Equation]^ that /<'w(P) reaches a maximum if = 0, a = 1, and a 
minimum is reached if a -0, /3 = 1. If the eigenvalues (curvatures) are not equal, we 
can show that the principle curvature directions are perpendicular. Eor this, we need 
to show that 5 is a self-adjoint operator. Thus, the equation {S v,w) = {v,5w) holds. 
We show this by using the property (n,^} — 0 and derive this with xj: 


dn d^f 

dxj ’ dxi ’ dxidxj 


) = 0 . 


We demonstrate that 5 is a self-adjoint operator with the basis 


df, , dn 5 /, , d^f , , dn df ^ .^,df, df ^ 


dxi ’ dxj 


dx: ’ dx ; 


’ dxjdxj 


dxj ’ dxi 


dxj ’ dxi 


Now, we show that the eigenvectors (principle curvature directions) are perpendic¬ 
ular if the eigenvalues (curvatures) are different: 


Ti<Vi,V2) = <5vi,V2) = <Vi,5V2) = T2<Vi,V2). 


The equation is only true if Vi,V 2 are perpendicular (and Ai + A 2 holds). 


3 Discrete Differential Geometry 

This section adapts the continuous differential geometry to discrete differential ge¬ 
ometry, the area of polygonal meshes that approximate continuous geometries. The 
following notation is used in the remainder of this paper. Let M c be a trian¬ 
gulated surface mesh. The mesh consists of vertices i e V with associated positions 
Pi e edges E = {(ij) \ ij e V}, and triangles T = e E). We 

write n/ as the normalized normal vector at vertex i. If nothing else is mentioned, we 
refer to normal vectors at vertices. Furthermore, Nil) denotes the neighbors of i. So, 
for every j € Nip), ii,j) € E holds. Furthermore, if we use a triangle for calculation, 
we always use this notation: given a triangle A = ii,j,k) with vertices p/,P;,pA:, and 
the edges are defined as Ci = P/ -Pj,e 2 = Pj -pA:,e 3 = Pi -p;. 
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mesh 


Fig. 5 The Voronoi diagram of different settings. In|(a)|a Voronoi diagram of a set of points is 


determined. In |(h)| the Voronoi area is calculated. If one of the triangles is obtuse, the area leaves 
the triangle. 


3.1 Voronoi Area 

We need to introduce the term Voronoi area, as it is important for the determination 
of the curvature. So, given are points in a 2D space. Every point is spread out in 
equal speed. If two fronts collide, they stop to spread out further at this region. 
After all fronts stopped, every point lies in a region that is surrounded by a front. 
This region is called a Voronoi region. Formally, given distinct points x, e in the 
plane, the Voronoi region for the point Xk is dehned as the set of points V{xk) with 

y(xi) = {x€R^ : ||x-Xi||<||x-X;||, 

See Figure [5(a)] for an example of a Voronoi diagram. To obtain the Voronoi area of a 
vertex on a surface mesh, the Voronoi area of each incident triangle is accumulated. 
The Voronoi area calculation is based on the method by Meyer et al. llMl . In case 
of a non-obtuse triangle, the Voronoi area at p,- is determined by the perpendicular 
bisector of the edges incident to p, . The point of intersection, the midpoint of the 
incident edges and the point itself dehne the endpoints of the Voronoi area. The 
triangle area of the Voronoi region equals: 

■5^a(Pi) = ^(lleill^ ■ cot(e2,e3) + HeslP ■cot(ei,e2)). 

In case of an obtuse triangle, the Voronoi area is equal half of the triangle area if the 
angle at p; is obtuse. Otherwise it is a quarter of the triangle area, see Figure [5(b)| 


3.2 Discrete Curvature 

The calculation of the curvatures as well as the principle curvature directions are 
important for a number of feature line techniques. Several approaches exist to ap- 
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proximate the curvatures. Some methods try to fit an analytic surface (higher order 
polynomials) to the mesh and determine the curvatures analytically ElEI- Another 
approach estimates the normal curvature along edges first and then estimates the 
shape operator Il8l fTSl [3111^ l46l . Other approaches are based on the calculation of 
the shape operator 5 Clllolllol. We use the curvature estimation according to M- 
After S is determined on a triangle basis, it is adapted to vertices. We already defined 
that S V yields the change of the normal in the direction of v: 

S \ - Dyll. 


This property is used to assess S for each triangle. When applying S to the edge ei, 
it should result in n, - ny because of the change of the normals along the edge. We 
need a basis of the tangent space of the triangle: 

ei . 62 

Cl = ] — 62 = - —:. 

|eil |e2l 

Afterwards, we build the orthogonal normalized basis vectors XA,yA by: 


Xa :=ei. 


XAX(e2XXA) 

||xAX(e2XXA)ir 


( 8 ) 


Applying the aforementioned property of the shape operator to all edges accord¬ 
ing to the basis leads to the following equation system: 


/<ei,XA)'] 


\<ei,yA)j 


/<e2,XA)'] 


\<e2,yA)j 


((e3,XA)) 


\<e3,yA)J 



see Figure for an illustration. Here, 
we have three unknowns (the matrix 
entries of the symmetric matrix S = 
and six linear equations.Thus, a 
least square method can be applied to 
fit the shape operator to approximate 
curvature for each triangle. Next, we 
need to calculate S for each vertex of 
the mesh. As the triangle basis nor¬ 
mally differs from each vertex tangent 
space basis, we need to transform the 
shape operator according to the new co¬ 
ordinate system. First, we assume that 
the normal Ha of the face is equal to the 


<n,-nj,XA)'\ 

<n,-nj,yA)/ 

{nj-nk,x^)\ 

<ni:-n;,XA)\ 
<ni--n;,yA)/ 



Fig. 6 The shape operator estimation is based on 
a local coordinate system, the edges and the nor¬ 
mals. 
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incident vertex normal n,. Hence, the basis (xajYa) of the triangle is coplanar to the 
basis (x,',y,) of the incident vertex i. Assuming we have the shape operator given in 
the vertex basis, then the entries can be determined by: 



As we have determined the shape operator in the basis (xA,yA), we can express 
the basis of the vertex by expressing the new coordinate system with the old basis 

X,- = axA +jSyA: 


Q' = <X,',Xa) 

jS = <x,',yA>- 


The entry ep can be determined by: 


ep^{a/3)s{^^. (10) 

The other entries can be calculated by analogous calculations. For the second case, 
we rotate the coordinate system of the triangle around the cross product of the nor¬ 
mal such that the basis of the vertex and the triangle are coplanar. Finally, we use this 
to determine the shape operator of the vertices. We determine the shape operators for 
all incident triangles of a vertex. Afterwards, we rotate the coordinate systems of the 
triangles to be coplanar with the basis of the vertex. Next, we re-express the shape 
operator in terms of the basis of the vertex. Then, we weight the shape operator ac¬ 
cording to the Voronoi area of the triangle and accumulate this tensor. Finally, we 
divide the accumulated shape operator by the sum of the weights. The eigenvalues 
provide the principle curvatures, and the eigenvectors give the principle curvature 
directions according to the basis. The pseudo-code [T] summarizes the algorithm. 

Please note that this algorithm can be generalized to obtain higher-order deriva¬ 
tives. It can be used to determine the derivative of the curvature as it is important for 
a specific feature line method. Formally, the derivative of the shape operator has the 
form: 


C = [DyS 




( 11 ) 
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Algorithm 1 Pseudo-code for the curvature estimation. 


for each triangle : 

Build basis accord, to Eq. [^ 

Determine S accord, to Eq.M 

for each vertex incident to the triangle: 
Rotate the triangle basis to the vertex basis 
Determine S in the new basis accord, to Eq.|10| 

Add this tensor weighted by the voronoi area 

end 

end 

for each vertex : 

Divide S by the sum of the weights 
Determine the eigenvalues and eigenvectors 

end 


For the determination of the change of the curvature in direction u, the tensor C has 
to be multiplied multiple times: 

DuA' = (u,(Z)v5 u Z)w5 u^ u). (12) 


3.3 Discrete Covariant Derivative 


First, we consider a linear 2D scalar field ^(x) - a-x\+li-X 2 +j and its gradient: 


SJ(f - 


\ikvj 



(13) 


To determine the gradient of a triangle A = {i,j,k) with scalar values (fi ^(p,), 
ipj iptpj), and ipk i^(Pi), we build a basis according to Equation[^and transform 
the points Pi,P;,Pjt 6 R.^ to p',p',p^ e by: 


P, = 



^/'<Pf-PLXA)\ 

\<Pf-PLyA>/ 


^/<P-i-PlXa)\ 

\<PA:-PLyA>/' 


This transformation describes an isometric and conformal map. The next step is 
a linearization of the scalar values ipi,ipj,ipk. We want to determine a scalar held 
ip' {n') - a ■ x'^ +I3- x' 2 +y such that 

(,tj'(p;) = ipi (,tj'(p') = ipj (^'(pp = ipk 

holds. These conditions yield the following equation system: 

[a p) (p' p} p^) + (r r r) = {v>i n) ■ 
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With p' = (q) we obtain the following solution; 


r = 

(a p) = [ipj - ifi ifik - V>}) (P} P^) ^ 


According to Equation the gradient of ip' is determined by 

The basis XA,yA yields the gradient in 

3D; 


V(^ = Q'-XA+j8-yA- 

Figure illustrates the gradient of a 
triangle. To determine the gradient per 
vertex, we use the same procedure as 
for the shape operator estimation. We 
transform the basis and weight the tri¬ 
angle gradient according to its propor¬ 
tion of the Voronoi area. 



3.4 Discrete Laplace-Beltrami Operator 

Several methods exist to discretize the Laplace-Beltrami operator on surface meshes. 
For an overview, we recommend the state of the art report by Sorkine ll43l . The op¬ 
erator can be presented by the generalized formula: 

^'/'(P/) = ^ Wij ((/:(P;) - ^(P,)). 
j 

Different weights wj give different discrete Laplace-Beltrami operators. For present¬ 
ing different versions of this operator it is preferable that it fulfills some properties 
motivated by the smooth Laplace-Beltrami operator; 

(Sym) The weights should be symmetric Wij — wji. 

(Loc) If (i,j) i E then Wjj - 0. 

(Pos) All weights should be non-negative. 

(Lin) If p; is contained in a plane and ip is linear, then Aipipi) - 0 should hold. 

In the following, we introduce different discrete Laplace-Beltrami operators. 

Combinatorial: For the combinatorial Laplace-Beltrami operator we have: 

fl, if(f,;)e£ 

Wij -1 

10 , otherwise. 
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Fig. 8 This figure illustrates the triangles with the angles for the weight calculation. 

This version may result in non-zero values on planar surfaces for linear scalar fields. 
Therefore, it violates (Lin). 

Uniform: Taubin 1471 suggested the uniform Laplace-Beltrami operator. The weights 
are determined by the number of neighbors of p,-; 

10 , otherwise. 

These weights also violate (Lin). 

Floater’s mean value: Floater iflTll proposed the mean value weights by the tangent 
of the corresponding angles: 



See Figure|^for the angles. These weights violate (Sym). 

Cotangent weights: MacNeal l2^ suggested the cotangent weights: 



cot(aij) + coiijiji), if (/, j) e E 
0 , otherwise. 


See Figure|^for the angles. On general meshes the weights can violate (Pos). 
Belkin weights: Belkin 13 suggested to determine weights over the whole surface: 



where A(a^) denotes the area of the triangle Ajt and h corresponds intuitively to the 
size of the neighborhood. This violates the (Loc) property. 






16 


Kai Lawonn and Bernhard Preim 



Fig. 9 In (a) the position of the zero crossing is determined and the points are connected. In 


(b) 


the isoline through a mesh is depicted. 


Results: The discussion leads to the question if there is any discrete Laplace- 
Beltrami operator which fulfills all required properties for an arbitrary surface mesh. 
Wardetzky et al. ifSOl showed that there is no such operator. The proof is based on a 
Schonhardt polytope which demonstrate that there is no Laplace-Beltrami operator, 
which does not violate any condition. 


3.5 Isolines on Discrete Surfaces 

For feature line methods, it is essential not to restrict the lines to the edges, as it 
is not desirable to perceive the mesh edges. Given is a surface mesh and a scalar 
field, we want to depict the zero crossing of the scalar field. Therefore, we linearize 
the scalar values for each triangle according to the values of the incident points. 
Afterwards, we look for points on an edge such that the linearized values of the 
scalar values of the connecting points are equal to zero. Having two points on two 
edges of a triangle, we connect them. Suppose we have a triangle with scalar values 
ipi > 0, ipj > 0 and ipk < 0. Thus, we know that somewhere on edge 62 and 63 there is 
a zero crossing. We determine t — and multiply t with edge 62 - This yields the 
position of the zero crossing on the edge. The position on the edge 63 is determined 
as well. Afterwards, both points will be connected, see Figure[9] 
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4 General Requirements of Feature Lines 

The generation of feature lines leads to several requirements, which have to be con¬ 
sidered for acquiring appropriate results. 

Smoothing: Most of the feature line methods use higher order derivatives. There¬ 
fore, the methods assume sufficiently smooth input data. For data acquired with 
laser scanners or industrial measurement process, smoothness cannot be expected. 
Discontinuities represent high frequencies in the surface mesh and lead to the gen¬ 
eration of distracting (and erroneous) lines. Several algorithms exist, which smooth 
the surface by keeping relevant features. Depending on the feature line method, dif¬ 
ferent smoothing algorithms can be applied. If the algorithm only uses the surface 
normals and the view direction, it is sufficient to simply smooth the surface normals. 
Geometry-based approaches, however, require to smooth the mesh completely. Op¬ 
erating only on scalar values, an algorithm which smoothes the scalar field around 
a certain region may be applied, too. 

Frame coherence: The application of feature line approaches or in general for non- 
photorealistic rendering makes it crucial to provide methods that are frame-coherent. 
This means, during the interaction the user should not be distracted by features that 
pop out or disappear suddenly. A consistent and continuous depiction of features 
should be provided in consecutive frames of animation. 

Filtering: Feature line algorithms may generate lines on salient regions as well as 
lines that result from small local irregularities, which may not be necessary to con¬ 
vey the surface shape or even annoying and distracting. Filtering of feature lines to 
set apart relevant lines from distracting ones is a crucial part of a feature line gener¬ 
ation. User-defined thresholds may control the rate of tolerance for line generation. 
Some algorithms use an underlying scalar field for thresholding. Lines are only 
drawn if the corresponding scalar value exceeds the user-defined threshold. Other 
methods integrate along a feature line, determine the value, and decide to draw the 
whole line instead of filtering some parts. We will also mention the filtering method 
of each presented feature line generation method. 


5 Feature Lines 

Line drawings were used extensively for medical visualization tasks, such as dis¬ 
playing tissue boundaries in volume data EHH, vascular structures 091 . neck 
anatomy ll25l and brain data ll^l45l . Furthermore, some higher order feature lines 
were qualitatively evaluated on medical surface data ca. The importance of fea¬ 
ture lines in medical visualization is discussed in OSlI . Feature line methods can be 
divided into image-based and object-based methods. Image-based approaches are 
not in the focus of this survey. These methods are based on an image as input. All 
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calculations are performed on the image with the pixels containing, for instance, an 
RGB or grey value. Usually, the image is convolved with different kernels to obtain 
the feature lines. The resulting feature lines are represented by pixels in the image 
space. These lines are mostly not frame-coherent. Comprehensive overviews of dif¬ 
ferent feature line methods in image space are given by This section 

presents selected object-based feature line methods. We will explain the methods 
and limitations. Further information on line drawings can be found in 1^1411 . 


5.1 Contours 


We refer to a silhouette as a depiction of the outline of an object as this is the 
original definition by Etienne de Silhouette. The contour is defined as the loci of 
points where the normal vector and the view vector are mutually perpendicular: 


<n,v) = 0, 



brain model 


where n is the normal vector and v is the view vector 
which points towards the camera. For the discrete case, 
we highlight edges as a contour whenever the sign of the 
dot product of the view vector with the normals of the 
incident triangle normals changes. The contour yields a 
first impression of the surface mesh. On the other hand, 
it is not sufficient to depict the surface well. The con¬ 
tour is not appropriate to gain a spatial impression of 
the object. Furthermore, it cannot depict salient regions, 
for instance strong edges. 

Summary: In the first place, the contour is necessary for gaining a first impression 
on the shape of the object. Unfortunately, spatial cues, as for instance strong edges, 
are not depicted. 


Fig. 10 The 
with contours. 


5.2 Crease Lines 

Crease lines are a set of edges where incident triangles change strongly. The di¬ 
hedral angle, i.e., the angle of the normals of the corresponding incident triangles, 
along the edges is calculated. The edge belongs to a crease line if the dihedral angle 
exceeds a user-defined threshold t. As the change of the normals is an indicator of 
the magnitude of the curvature, one can state that all points contribute to a feature 
line if the underlying absolute value of the maximum curvature exceeds a threshold: 


Ki>T or (n,',nj)>T', 
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for adjacent triangles with corresponding normals ni,nj. 
Afterwards, all adjacent vertices which fulhll the prop¬ 
erty are connected. These feature lines need to be com¬ 
puted only once, since they are not view-dependent. Fur¬ 
thermore, these lines are only drawn along edges. 
Summary: Crease lines display edges where the dihe¬ 
dral angle is large. Strong edges are appropriately de¬ 
picted, but if the object has small features, this method 
is not able to depict only important edges. This is caused 
by the local determination of the dihedral angle without 
concerning a neighborhood. Even smoothing the surface 
mesh would not deliver proper line drawings. Further¬ 
more, this method is only able to detect features on edges 



Fig. 11 The brain model 
with crease lines and con¬ 
tours. 


5.3 Ridges and Valleys 


Ridges and valleys were proposed by Interrante et 
ai. im and adapted to triangulated surface meshes by 
Ohtake et al. llT5l . These feature lines are curvature- 
based and not view-dependent. The computation is 
based on the principle curvature ki as well as the as¬ 
sociated principle curvature direction ki with |a:i| > \k 2 \. 
Formally, ridges and valleys are dehned as the loci of 
points at which the principle curvatures assume an ex¬ 
tremum in the principle direction: 


DkiKi - 0 . 



Fig. 12 The brain model 
with ridges and valleys, and 
contours. 


According to two constraints, the sets of points are 
called 


DkiDk^Ki 


I < 0, and Ki > 0: ridges 
> 0, and Ki < 0: valleys. 


(14) 


To determine the ridge and valley lines, we hrst need to compute the principle cur¬ 
vatures and their associated principle curvature directions, recall Section[T^ After¬ 


wards, we determine the gradient of ki for each vertex, see Section 3.3 Finally, we 


compute the dot product of the gradient and the associated principle curvature direc¬ 
tion ki. This yields the scalar value of for each vertex. Next, we distinguish 

between ridges and valleys and determine Dt, Dk^ ki for each vertex. Here, we need 
again the gradient of each vertex with the value Dtj ki and determine the dot prod¬ 
uct of the result with ki. Hence, we gain two scalar values per vertex: D^iKi and 
DkiDk^Ki. Afterwards, we assess the zero-crossing of the hrst scalar value, recall 
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Section [X5| We connect the zero crossings in every triangle for which one condition 
of Equationholds. The filtering of the lines is again performed by employing an 
user-defined threshold. The integral along each ridge and valley line is determined 
according to the underlying curvature. If the magnitude of the integral exceeds the 
threshold for ridges or valleys, the line is drawn. 

Summary: The calculation is solely based on the curvature and therefore view- 
independent. This method is able to detect small features. The filtering depends on 
the underlying curvature and the length of the curve. Therefore, a long line with 
small curvature has also the chance to be drawn as a small line with high curvature. 
This strategy emphasizes also long feature lines. Ridges and valley lines are very 
susceptible to noise, since this method is of 3rd order. Therefore, small discontinu¬ 
ities on the surface mesh lead to erroneous derivatives and this error propagates for 
each further derivative. A crucial task for this method is to guarantee a smoothed 
mesh to obtain reasonable results. From an artist’s point of view, some features may 
be more highlighted than others from different points of view. This is caused by the 
different perception of an object and by various light positions. For this task, the 
ridge and valley lines are not appropriate due to the restriction of view-independent 
results. 


5.4 Suggestive Contours 


Suggestive contours are view-dependent feature lines in¬ 
troduced by DeCarlo et al. IfTTlI . They extend the defini¬ 
tion of the contour. These lines are defined as the set of 
minima of (n, v) in the direction of w, where n is the sur¬ 
face normal, v is the view vector which points towards 
the camera, and w = (Id - nn^)v is the projection of the 
view vector on the tangent plane. Formally: 


D„(ii,v) = 0 and D„D„(n,v)>0. 



Fig. 13 The brain model 

Another equivalent definition of the suggestive contours suggestive contours and 

is given by the radial curvature Kr- It is defined as the contours, 
curvature in direction of w. As seen in Equation]^ this 
curvature can be determined by knowing the principle 

curvature directions as well as the corresponding curvatures. Therefore, the defini¬ 
tion of the suggestive contours is equivalent to the set of points at which the radial 
curvature Kr is equal 0 and the directional derivative of Kr in direction w is positive: 


Kr -0 and Dy,Kr > 0. 


The filtering strategy is to apply a small threshold to eliminate suggestive contour 
points where the radial curvature in direction of the projected view vector is very 
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low. Additionally, a hysteresis threshold is applied to increase granularity. 
Summary: Suggestive contours extend the normal definition of the contour. This 
method depicts zero crossing of the diffuse light in view direction. This can be seen 
as inflection points on the surface. This method is of 2nd order only and thus less 
susceptible to noise. Unfortunately, suggestive contours are not able to depict some 
sorts of sharp edges, which are in fact noticeable features. For instance, a rounded 
cube has no suggestive contours. 


5.5 Apparent Ridges 


Apparent ridges were proposed by Judd et al. E^ . 
These feature lines extend the definition of ridges by 
a view-dependent curvature term. Therefore, a projec¬ 
tion operator P is used to map the vertices on a screen 
plane V. The orthonormal basis of the screen plane is 
given by (vi,V 2 ). Assume we have a parametrized sur¬ 
face /: / c —> R^. Then the projection of / onto V is 
given by: 


F(x) = 


/<vi,/(x))\ 

l<V2,/(x))j- 



Fig. 14 The brain model 
with apparent ridges. 


The Jacobian Jp of P can be expressed as: 




In the discrete case with surface meshes, the Jacobian can be expressed by a basis 
for the tangent plane ( 61 , 62 ): 


j ^/<vi,ei) (¥1,62)) 

^ \<V2,ei) <V 2 , 62 )/' 

If a point p' on the screen plane is not a contour point, there exists a small neigh¬ 
borhood where the inverse of P exists. Normal vectors n' at a point p' on the screen 
plane are defined as n'(p') := n(P“*(p')). The main idea is to build a view-dependent 
shape operator 5' at a point p' on the screen as 


5'(w') = £>w'n' 


where w' is a vector in the screen plane. The view-dependent shape operator is 
therefore defined as: 


S' ^SJp\ 
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Here, the basis of the tangent space expressing S and Jp must be the same. In con¬ 
trast to the shape operator, the view-dependent shape operator is not a self-adjoint 
operator, recall Section [ZS] Therefore, it is not guaranteed that 5' has two eigenval¬ 
ues, but it has a maximum singular value /cj: 

K. - max ||5'(w')ll- 

||w|| = l 


This is equivalent to hnd the maximum eigenvalue of S'^S' and to take the 
square root. The corresponding singular eigenvector t' is called the maximum view- 
dependent principle direction. The rest of the method is similar to the ridge and 
valley methods. Formally, apparent ridges are defined as the loci of points at which 
the view-dependent principle curvature assumes an extremum in the view-dependent 
principle direction; 


Dt'K[ =0 and DpDi>K[ < 0. 


The sign of k' is always positive. To distinguish between ridge lines and valley lines, 
we may compare the sign of the object-space curvature: 

I < 0, ridges 
Kl j 

1 >0, valleys. 

The calculation of the directional derivative is different from the other methods. This 
calculation is performed with hnite differences. Therefore, we transform the singu¬ 
lar eigenvector t' to object space t using the corresponding basis of the associated 
vertex i. Furthermore, we need the opposite edges of the vertex and determine two 
points wi, W 2 on the edges such that t and the edges are orthogonal and wi, W 2 are 
the dropped perpendiculars of t to the corresponding edges. The directional deriva¬ 
tives are determined by averaging the hnite differences of the curvatures between 
Pi and wi, W 2 . The curvature of wi, W 2 is assessed by linear interpolation of the 
endpoints of the associated edge. Having the principle view-dependent curvature 
direction t', we need to make it consistent over the mesh because it is not well- 
dehned. Therefore, t' is hipped in opposite direction whenever it does not point to 
the direction where the view-dependent curvature is increasing. The zero-crossings 
are determined by checking if the principle view-dependent curvature directions of 
the vertices along an edge point are in the same direction. Only in this case there is 
no zero-crossing. Pointing in different directions means that the enclosing angle is 
greater than 90 degrees. The zero crossing is determined by interpolating the values 
of the derivatives. To locate only maxima, a perpendicular is dropped from each ver¬ 
tex to the zero crossing line. If the perpendiculars of the vertices of an edge make 
an acute angle with their principle view-dependent curvature directions, the zero 
crossing is a maximum. Otherwise, the zero crossing is a minimum. To eliminate 
unimportant lines, a threshold based on the view-dependent curvature is used. 
Summary: Apparent ridges incorporate the advantages of the ridges and valley lines 
as well as the view dependency. They extend the ridge and valley dehnition by in- 
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troducing view-dependent curvatures. This method is able to depict salient regions 
as sharp edges. Unfortunately, the 3rd order computation leads to low frame rates 
and to visual clutter if the surface mesh is not sufficiently smoothed. 


5.6 Photic Extremum Lines 


Photic extremum lines (PELs) were introduced by Xi et 
al. II 52 I . These feature lines depict regions of the surface 
mesh with significant variations of illuminations. This 
method is based on the magnitude of the light gradi¬ 
ent. Formally, these lines are defined as the set of points 
where the variation of illumination along its gradient di¬ 
rection reaches a local maximum: 


Z)„||V/|| = 0 and DMl^fW < 0, 



with w = jj^. Normally, / is used as the headlight il¬ 
lumination: / := <n, v) with n as the normal vector and 
V as the view-vector. PELs have more degrees of free¬ 
dom to influence the result by adding more light sources. Thus, the scalar value of 
/ changes by adding the light values of the vertices by other lights. Noisy photic 
extremum lines are filtered by a threshold which is based on the integral of single 
connected lines. The strength T of a line with points xq, ... ,x„ is determined by: 


Fig. 15 The brain model 
with photic extremum lines. 


IIV/II = Xi 

1=0 


I|V/(X;)|| + ||V/(X,-+1)|I 


l|x,-x,+i| 


If T is less than a user-defined threshold, the line is canceled out. 

Summary: Photic extremum lines are strongly inspired by edge detection in image 
processing and by human perception of a change in luminance. It uses the variation 
of illumination. The result may be improved by adding lights. Beside the filtering 
strategy to integrate over the lines and accumulate the magnitude of the gradient, 
the noise can also be reduced by adding a spotlight that directs to certain regions. 
Nevertheless, smoothing is necessary to gain reasonable results. Here, the smooth¬ 
ing of the normal is sufficient as the computation is mainly based on the normals. 
However, the computation has high performance costs. The original work was im¬ 
proved by Zhang et al. Il54l to significantly increase the runtime. 




24 


Kai Lawonn and Bernhard Preim 


5.7 Demarcating Curves 


Demarcating curves were proposed by Kolomenkin et 
ai. m. These feature lines are defined as the transition 
of a ridge to a valley line. To determine these lines, the 
derivative of the shape operator has to be calculated, re¬ 
call Equation [TT] 

C = (Z)v5 D^S). 

The demarcating curves are defined as the set of points 
where the curvature derivative is maximal; 


{w,Sw} — 0 with w = arg max Dy/r. 

I|V||=1 



Fig. 16 The brain model 
with demarcating curves and 
contours. 


The values for w can be analytically found as the roots 

of a third order polynom. This is obtained by setting v = (co"(e)) combining this 
with Equation [T^ A user-defined threshold eliminates demarcating curves, if it ex¬ 
ceeds the value of D^k. 

Summary: Demarcating curves are view-independent feature lines displaying re¬ 
gions where the change of the curvature is maximal. Therefore, higher-order deriva¬ 
tives are used. A 2x2x2 rank-3 tensor is determined. This method can be used 
to illustrate bumps by surrounding curves. The advantage of the method is to en¬ 
hance small features. Especially when combined with shading, this approach has 
its strength in illustrating archaeology objects where specific details are important, 
e.g., old scripts. For this application, view-dependent illustration techniques are not 
recommended because details need to be displayed for every camera position. Con¬ 
trary, due to higher-order derivatives, the method is sensitive to noise and is not well 
suited for illustrative visualization. 


5.8 Laplacian Lines 

Laplacian lines were proposed by Zhang et al. Il55l . The introduction of these lines 
was inspired by the Laplacian-of-Gaussian (LoG) edge detector in image process¬ 
ing and aims at a similar effect for surface meshes. The idea of the LoG method 
is to determine the Laplacian of the Gaussian function and to use this kernel as a 
convolution kernel for the image. Laplacian lines calculate the Laplacian of an il¬ 
lumination function / and determine the zero crossing as feature lines. To remove 
noisy lines, the lines are only drawn if the magnitude of the illumination gradient 
exceeds a user-defined threshold r; 


zl/ = 0 and IIV/II > T, 




Feature Lines for Illustrating Medical Surface Models 


25 


where A is the discrete Laplace-Beltrami operator on the surface mesh and / is the 
illumination with / := (n,v). Here, the discrete Laplace-Beltrami operator with the 
Belkin weights is used, as introduced in Section 3.4 The advantage of this method 
is the simplified representation of the Laplacian of the illumination; 


/f/(p) =zl<n,v) 

= <^n,v). 

Here, An is the vector Laplace operator in the Euclidean space. 


This is just a composite of the Laplacian of the differ¬ 
ent components. Thus, the algorithm consists of a pre¬ 
processing step to calculate the Laplace-Beltrami oper¬ 
ator with the Belkin weights of the components of the 
normal An. During runtime, the algorithm detects the 
zero crossings of (An,v) and checks if the magnitude of 
||V/|| exceeds the user-defined threshold. 

Summary: The Laplacian lines are strongly inspired 
by edge detection algorithms in image processing. This 
method is based on the Laplacian-of-Gaussian. Basi¬ 
cally, the method searches for zero crossings in the 
Laplacian of the illumination. The computational effort 
can be simplified by a preprocessing step. Thus, interactive frame rates for geo¬ 
metric models of moderate size are possible during the interaction. Similar to other 
higher order methods, this approach also assumes well smoothed surface normals. 
The Belkin weights for the Laplace-Beltrami operator have a smoothing effect for 
the Laplacian line generation. This method illustrates sharp edges well, but is not 
suitable for round corners. 



Fig. 17 The brain model 
with Laplacian lines. 


6 Discussion and Comparison 

This section deals with general proper¬ 
ties of the different feature line meth¬ 
ods. We discuss the different approaches 
to derive first recommendations which 
method may be used for which kind 
of geometry. First, we list all feature 
line methods in Table [T] and name 
different properties and the order of 
the corresponding method. Furthermore, 
in Figure the higher-order feature 
lines are illustrated on an analytic func¬ 
tion. 


Name 

Order 

View-dep. 

Contours 

1 

yes 

Crease Lines 

1 

no 

Ridges & Valleys 

3 

no 

Suggestive Contours 

2 

yes 

Apparent Ridges 

3 

yes 

Photic Extremum Lines 

3 

yes 

Demarcating Curves 

3 

no 

Laplacian Lines 

3 

yes 


Table 1 List of different feature line meth¬ 
ods with derivative order and view depen¬ 
dency. 
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A f- A 1» TT 

/\ 

(a) Ridges and Valleys, Ap- (b) Suggestive Contours, (c) Photic Extremum Lines, 
parent Ridges Demarcating Curves Laplacian Lines 


Fig. 18 Drawing of an analytic function with illustrated feature line positions. In (a) the ridges 


are denoted in orange and the valleys are illustrated in cyan. For this function with fixed view 
direction, the apparent ridges coincide with ridge and valley lines. In |(b)| the suggestive contours 
and the demarcating curves are the same. In |(c)| the photic extremum lines and the Laplacian lines 
coincide. 


The benefit of feature lines is motivated by the visual perception. In l30l it is stated 


Name 

Sharp Edges 

Round Edges 

Bumps (s.w.) 

Bumps (top) 

Contour 

Deformation 

Contour 

El 

E 

E 

E 

Ef 

Ef 

Crease Lines 

[Zf 

E 

E 

E 

E 

Ef 

Ridges & Valleys 

[Zf 

tzf 

E 

Ef 

E 

E 

Suggestive Contours 

El 

E 

Ef 

Ef 

E 

Ef 

Apparent Ridges 

EZf 

tzf 

eT 

Ef 

Ef 

E 

Photic Extremum Lines 

[Zf 

E 

E 

Ef 

Ef 

Ef 

Demarcating Curves 

El 

E 

E 

Ef 

E 

E 

Laplacian Lines 

Ef 

E 

E 

Ef 

Ef 

E 


Table 2 List of supported feature by the methods. The different features are illustrated in Figure|l9| 


that the first stage of the assessment of the shape is done by extracting features, such 
as contours. These characteristics help to understand the shape. The illustration of 
shapes with feature lines cannot be seen as an alternative to shading. It is rather 
an additional concept. Kolomenkin et al. 1241 showed that their demarcating lines 
support the shading and can extract text from archaeology objects. However, for ex¬ 
amining structures where the whole object inherits important information, feature 
lines should not be used solely. For data where the scene can be divided into focus 
and context objects, feature lines can be applied to the context objects. Furthermore, 
feature lines can also be used to enhance focus with additional shading. 

Depending on the underlying model, we may recommend different techniques. 
Most of the feature lines are able to depict the contour, but this depends strongly 
on the bending of the surface at the contour. Especially apparent ridges and photic 
extremum lines are able to draw contour lines, but in our experiments we noticed 
that activating the contour enhances the visual impression because some parts of the 
contour were missing. If the surface model is an assembly with sharp edges we rec¬ 
ommend to use ridge and valley lines or apparent ridges. These features often appear 
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in medical models models like implants or prostheses. For simple models with only 
a few sharp edges, crease lines may be appropriate as well. If the models have a lot 
of round edges the answer for the right feature line method is a matter of taste. These 
features appear in models like vascular surfaces or organs. For scenarios where it is 
important to illustrate details for every camera position, we recommend ridges and 
valleys as well as demarcating curves. From an artistic point of view, suggestive 
contours, photic extremum lines, and Laplacian lines should be chosen. The reason 
for this suggestion is that especially for a rounded cube the photic extremum lines 
and Laplacian lines generate double lines around the feature to denote the rounded 
edge. If the user wants to visualize the line along the edge, the ridge and valley lines 
or the apparent ridges should be used. For this case, the crease line approach is not 
useful because it depicts only edges with specific greater value of the dihedral angle. 
Therefore, too many lines may be generated. If the surface has many crevices, we 
recommend the suggestive contours. They illustrate the inflection points of valleys. 
Round corners are often represented in many organic structures like livers or bones, 
see Figurel^for a femur model or a skull model. 

Table |2]lists all possible features and Figure shows the different features. 
Please note that the assessment of the suitability of a method - marked in the table 
- necessarily is a subjective assessment by the authors and two artists. For instance, 
regarding the property whether the methods are able to detect round edges, we mean 
if it detects the specific round feature. As already mentioned, it does not reflect the 
ability to enhance the round edge from an illustrative or artist point of view. This 
concerns the ability to depict bumps. In agreement with artist, the bump shown from 
a sideway (s.w.) perspective would be illustrated such that it depicts the smooth 
transition from the ground to the dent. The drawing of the surrounding circle of 
the bump is not desirable as it conveys a sharp transition from the bump to the 
ground. For bumps shown from the top perspective it is sufficient if a round circle 
is drawn. Bumps can occur as polyps or blebs on a cerebral aneurysm. Especially 
blebs are important anatomical features to be detected because they are an indicator 
for rupture. Blebs can also occur as polyps in CT colonography. 

We also listed the property deformation in the table. This characteristic means if 
the corresponding method is able to illustrate the features of deformable surfaces, 
e.g., animated objects, in real-time. As an example, Oeltze et al. 041 analyzed my¬ 
ocardial perfusion data. The focus lies on the examination of the infarction scarf 
on the left ventricle. In this paper, the left ventricle is illustrated as context infor¬ 
mation. Using the time-dependent data, it would also be possible to illustrate the 
context information with some feature line methods during the animation. 

For example suggestive contours have two definitions of how to assess the feature 
lines. One is curvature-based and the other is light-based. With the second defini¬ 
tion, no preprocessing is needed to assess the curvature and the principle curvature 
directions. This is in contrast to ridges and valleys and apparent ridges. Therefore, 
these algorithms are not able to compute the feature lines during the deformation. 
Photic extremum lines are also able to compute the feature lines during runtime be¬ 
cause of the light and view dependency. The Laplacian lines need to precompute the 
Laplacian of the normals. Hence, this method is not suited for deformations. 
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Fig. 19 Different surface features are illustrated with shading (SH) and in different high-order 
feature line methods: suggestive contours (SC), ridge and valley (RV), apparent ridges (AR), photic 
extremum lines (PEL), demarcating curves (DC), and Laplacian lines (LL). 


Furthermore, Figure 20 shows some exemplary models illustrated with higher 
order feature lines. Three typical models in the discrete differential geometry field 
(cow, Buddha, Max Planck) as well as three models from the medical image data 
(brain, femur, skull) are presented. 

In summary, current feature lines are not suitable for the depiction of anatomical 
structures directly derived from medical image data because the underlying surfaces 
are too noisy. Advanced smoothing algorithms are necessary to reduce artifacts, but 
preserve important anatomical structures. For the depiction of a sparse representa¬ 
tion of the model in a context-aware manner, the feature line methods can be used. 


6.1 Medical Application 

As stated, feature line methods can be seen as an illustrative visualization method 
that can enhance shading or as an alternative in the focus-and-context visualization. 
In this section, we list different application fields where illustrative visualization 
is useful for effectively depicting medical data. At the end, we list possible helds 
where feature lines can be used to encode context information. 

Fischer et al. na proposed to use illustrative visualizing tools to depict structures 
of hidden surfaces. The rendering style is tailored for understanding spatial relation¬ 
ships and for visualizing hidden objects. Born et al. 0 used illustrative techniques 
to depict stream surfaces. Their techniques are very useful for the visualizing of 
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SH SC RV AR PEL DC LL 


Fig. 20 Selected models depicted in shading and higher-order feature lines. 



(a) (b) (c) 




Fig. 21 Different medical application fields where feature lines can be used to illustrate surround¬ 
ing objects. 


complex flow structures. In the area of brain data, Jainek et al. Il22ll suggested to 
use a hybrid visualization method to illustrate mesh and volume rendering. Their 
approach is efficient for the exploration in clinical research. Chu et al. 13 proposed 











30 


Kai Lawonn and Bernhard Preim 


a guideline of various rendering techniques. They combined, e.g., isophote-based 
line hatching and silhouette drawing, for illustrative vascular visualization. Different 
rendering techniques for medical applications were presented by Tietjen et al. EH. 
An example for illustrative visualization for liver surgery can be found in lfT9ll . 

GlaBer et al. ifThI presented an approach to visualize 3D cluster results. Here, 
the medical researcher can analyze the whole 3D scene with different cluster results 
where he can also select interesting objects. The surrounding objects become con¬ 
text information. Thus, we propose to illustrate them with feature lines. In this case. 


illustrates the main object with unselected objects illustrated with feature lines. 

In the field of endoscopic views, the identification of polyps is necessary. Once 
the polyps are detected, they can be illustrated in such a way that the endoscopic 
views are used for context information. In Figure [21(b)[ we used suggestive contours 
for the vessel and diffuse shading for the polyps. 

In Figure [2l(c)| we visualized the portal vein and three liver segments. The portal 
vein is illustrated in diffuse shading in red. The liver segments are visualized in 
diffuse shading with transparency and photic extremum lines. 


we used the contour because the objects does not inherit much features. Figure 21(a) 


7 Conclusion 

We have summarized the most common feature line methods for object space-based 
presentations of 3D meshes as they are frequently used in medicine and molecular 
biology. The presentation of the different methods was also covered by two basic 
sections. We did not only list the most common feature line methods and their cal¬ 
culation in the discrete space, but also provided the mathematical background to 
explain the calculation from the differential geometry point of view. Our goal was 
to present an extensive list of feature lines on the one hand, and to equip the reader 
with basic knowledge of differential geometry on the other hand. The graduated 
student may be able to follow the different methods and to implement every feature 
line algorithm based on our explanations in the field of discrete differential geom¬ 
etry. Therefore, our survey and tutorial may also be used by students who are new 
in the field of illustrative rendering. Furthermore, this survey may also be used as 
a starting point for the development of new feature line methods. The potential of 
advanced and recently introduced feature line techniques is currently not exploited 
in the display of medical surface models. The careful application of these methods 
and perception-based evaluations are left open for future work. 
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